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ABSTRACT 

The weak-lensing science of the Large Synoptic Survey Telescope (LSST) project drives the need 
to carefully model and separate the instrumental artifacts from the intrinsic shear signal caused by 
gravitational lensing. The dominant source of the systematics for all ground-based telescopes is the 
spatial correlation of the point spread function (PSF) modulated by both atmospheric turbulence and 
optical aberrations in the telescope and camera system. In this paper, we present a full field-of-view 
simulation of the LSST images by modeling both the atmosphere and the system optics with the most 
current data for the telescope and camera specifications and the environment. To simulate the effects 
of atmospheric turbulence, we generated six-layer Kolmogorov/von Karman phase screens with the 
parameters estimated from the on-site measurements. LSST will continuously sample the wavefront, 
correcting the optics alignment and focus. For the optics, we combined the ray-tracing tool ZEMAX 
and our simulated focal plane data to introduce realistic residual aberrations and focal plane height 
variations. Although this expected focal plane flatness deviation for LSST is small compared with that 
of other existing cameras, the fast f-ratio of the LSST optics makes this focal plane flatness variation 
and the resulting PSF discontinuities across the CCD boundaries significant challenges in our removal 
of the PSF-induced systematics. We resolve this complication by performing principal-component- 
analysis (PCA) CCD-by-CCD, and interpolating the basis functions derived from the analysis using 
conventional polynomials. We demonstrate that this PSF correction scheme reduces the residual PSF 
ellipticity correlation below 10 -7 over the cosmologically interesting (dark matter dominated) scale 
10' — 3°. From a null test using Hubble Space Telescope (HST) Ultra Deep Field (UDF) galaxy images 
without input shear, we verify that the amplitude of the galaxy ellipticity correlation function, after 
the PSF correction, is consistent with the shot noise set by the finite number of objects. We conclude 
that the current optical design and specification for the accuracy in the focal plane assembly are 
sufficient to enable the control of the PSF systematics required for weak-lensing science with LSST. 

Subject headings: Astronomical Instrumentation — gravitational lensing — dark matter — cosmology: 
observations — Astrophysical Data — Astronomical Techniques 



1. INTRODUCTION 

The Large Synoptic Survey Telescope (LSST) has been 
designed to provide a deep six-band (0.3-1.1 ^m) astro- 
nomical imaging survey of more than 20, 000 square de- 
grees of the southern sky. Using active optics, the 8.4- 
meter aperture and 9.6 square degree field-of-view of the 
telescope will allow ~ 1000 visits to each patch of sky 
in ten years with a final depth reaching r ~ 27.5 and a 
mean delivered PSF better than 0.7" (LSST Science Col- 
laborations et al. 2009). The LSST project will deliver 
fully-calibrated, science quality images, catalogs, and de- 
rived data products to the US public with no proprietary 
period. We refer readers to LSST Science Book II (LSST 
Science Collaboration 2009) for extensive discussions on 
the science that the survey will enable. 

Among the most critical science applications that the 
LSST will revolutionize is weak gravitational lensing. 
The key observable in weak gravitational lensing is the 
shape of distant galaxies weakly distorted by foreground 
masses. Because each galaxy is sheared only by a small 
amount, the lensing signal must be extracted from an 
average over a population of galaxies. Consequently, the 
total number of resolved galaxies at the usable surface 
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brightness limit is one useful measure to quantify the 
effectiveness of an astronomical survey for weak-lensing 
science. Another important factor is the number den- 
sity of usable galaxies (shear noise weighted effective 
number), which determines the lower limit of the an- 
gular scale that lensing can probe. The LSST project 
is designed to maximize the scientific return balancing 
both the total number (area) and the number density 
within the ten-year nominal mission period. A sky cov- 
erage of approximately 20, 000 square degrees is re- 
quired to keep the sample variance below the statisti- 
cal limit, and the high-resolution (< 0.7" median deliv- 
ered seeing), deep imaging (r ~ 27.5), providing > 40 
high S/N galaxies per square arcmin enabling a reli- 
able probe of matter structure on galaxy cluster scales. 
Residual uncorrected PSF shape variation requirements 
are summarized in the Science Requirement Documents: 
http://www.lsst.org/filesdocs/SRD.pdf] 

These performance requirements determine the desired 
level of accuracy in point spread function (PSF) cor- 
rection. The effect of the PSF on the galaxy shape 
measurement cannot be overemphasized. Anisotropic 
PSFs mimic gravitational lensing signal, generating a co- 
herent alignment of galaxy shapes. Even without this 
anisotropy, PSFs circularize shapes of barely resolved 
galaxies and dilute the intrinsic lensing signal. For 
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various reasons, the observed PSF changes in size and 
anisotropy across the field and also over time. Moreover, 
we need the PSF extrapolated to the positions of galax- 
ies, whereas the PSF model is constructed from the high 
S/N stars that are much less densely distributed. There- 
fore, in addition to very low surface brightness limits, a 
key factor in precision lensing lies in one's ability to care- 
fully model and remove the complicated effects of PSFs. 

Although the seeing at the LSST site is among the 
best in existing ground-based facilities, the accurate de- 
scription of the PSF is highly challenging because of 
the optical and focal plane design. As described below, 
the telescope optics and camera are continuously aligned 
via wavefront sensing. Optimized for the unprecedented 
large field of view, the effective /-ratio of LSST is small 
// ~ 1.2, which makes the optical aberration highly 
sensitive to alignment and focal plane deviation. Be- 
cause LSST's focal plane will be tiled with 189 4k x 4k 
CCDs, any height fluctuation both across and within the 
CCDs will translate into a very complicated PSF pattern, 
which is characterized by abrupt changes across the CCD 
gaps and smooth, but possibly high-frequency variations 
within the CCDs. Any sub-optimal modeling of these 
PSFs will leave systematic residuals on the scales of the 
CCD sizes and the "potato chip" effeclr] These system- 
atic residuals, mimicking lensing signals^ will, of course, 
hamper the correct interpretation of the lensing analysis. 

Unfortunately, with the existing algorithms the ability 
to precisely describe the PSF alone does not guarantee 
the success in the extraction of gravitational shear to the 
accuracy that future LSST-like weak-lensing surveys re- 
quire. This obstacle is well noted in the recent large 
collaborative shear measurement campaign GREAT08 
Challenge (Bridle et al. 2010). The campaign let the 
participants perform blind shear measurements on simu- 
lated images, where the PSF was known, but the input 
shear was unknown. Although many algorithms were 
shown to provide the accuracy for high S/N images suit- 
able for the existing survey data, no method came close 
to the target accuracy Q ~ 1000 when the noise level 
of the images matched the realistic value (see Bridle et 
al. 2010 for the definition of Q). Nevertheless, it is im- 
portant to note that the main challenge is purely math- 
ematical/statistical and thus improvable as more mathe- 
matical sophistication is incorporated in the algorithms. 
For example, Bernstein (2010) claims that when the bias 
arising when the true galaxy profiles do not match the 
models being fit is properly addressed, the resulting al- 
gorithm can achieve Q ~ 3000 for high S/N galaxies of 
the GREAT08 sample. 

We launched the LSST shear calibration project in or- 
der to diagnose the key factors in the telescope and cam- 
era engineering specifications affecting the weak-lensing 
science and to develop new algorithms for optimal shear 
extraction in the presence of different combinations of 
systematics. The current paper is the first in the series 
of this topic with an emphasis on the realization, char- 
acterization, and reconstruction of LSST PSFs. 

The paper is organized as follows. In fj2] we provide a 
brief review on the telescope optics. Our implementation 
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of the atmospheric turbulence is described in |3j In Sec- 
tion S|4j we discuss the focal plane design of LSST and 
the resulting behavior of the PSF. ^presents our algo- 
rithm to measure and reconstruct the PSFs. In f|6j we 
describe our simulation. Finally, the detailed comparison 
between the observed and modeled PSFs is presented in 
Sj7] before the conclusion in Sj8] 

2. LSST ENVIRONMENT AND OPTICS 

Figure [T| shows the optical design of the LSST, which is 
a modified Paul-Baker three-mirror system. The active 
optics telescope consists of an 8.4 m f/1.18 concave pri- 
mary, a 3.4 m f/1.0 convex secondary, and 5.2 m f/0.83 
concave tertiary mirrors. The final focal ratio f/1.23 (fo- 
cal length of 10.3 m) gives a plate scale 0.0197"//Ltm over 
the 64 cm diameter focal plane (3.5° x 3.5°). The 27% 
obscuration yields a total effective light collecting area 
of 35 m 2 , which corresponds to an effective 6.7 m di- 
ameter clear circular aperture. Current wide field tele- 
scopes were not designed for the stringent PSF demands 
of LSST weak lensing goals. 

This optical design provides a wide field of view while a 
maintaining uniform image quality across the field (Sep- 
pala 2002). The remarkably small variation of the PSF 
size over the field is illustrated in Figure [2j where the 
diffraction-limited PSF images are generated with ZE- 
MAX at the nominal flat focal plane. The peak-to-valley 
80% encircled energy radius is within ~ 7% of the mean 
value from the center to ~ 1.4° for simulated i filter point 
sources. The maximum deviation (~ 17%) is seen at the 
edge of the field (1.75°). The image brightness is stable 
across the field, providing a nearly uniform illumination 
within a radius of 1.2° and about 10% decrease at the 
edge (Seppala 2002; LSST Science Collaborations 2010). 
Finally, the LSST field distortion is remarkably small: 
less than 0.1% over the full field. 

LSST will have continuous correction of its optics, us- 
ing curvature wavefront sensing. Four special purpose 
rafts, mounted at the corners of the science array, contain 
wavefront sensors and guide sensors (right panel of Fig- 
ure fl]). Wavefront measurements are accomplished using 
curvature sensing, in which the spatial intensity distribu- 
tion of de-focused stars is measured at equal distances on 
either side of focus. Each curvature sensor is composed 
of two CCD detectors, with one positioned slightly above 
the focal plane, the other positioned slightly below the fo- 
cal plane. The CCD technology for the curvature sensors 
is identical to that used for the science detectors in the 
focal plane, except that the curvature sensor detectors 
are half-size so they can be mounted as an in-out defo- 
cus pair. Detailed analyses (Manuel et al. 2010) have 
verified that this configuration can reconstruct the wave- 
front within the < 0.2/xm accuracy. These four corner 
rafts also hold two guide sensors each. The guide sensors 
monitor the locations of bright stars at a frequency of 
~ 10 Hz to provide feedback for a loop that controls and 
maintains precision tracking of the telescope during an 
exposure. 

The fast focal ratio and the rapid pointing changes 
make any hardware-based atmospheric dispersion cor- 
rection technically difficult. Consequently, the effect of 
the atmospheric dispersion sets the maximum angle away 
from the zenith. We estimate the atmospheric differen- 
tial refraction for the six proposed filters (it, g, r, i, z, and 
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Y) using the models summarized in Filippenko (1982). 
Figure [3] displays the results for the input parameters of 
/ = 8 mm Hg (water vapor pressure), T — 5°C (atmo- 
sphere temperature), and P — 520 mm Hg (atmospheric 
pressure), which are the approximate average conditions 
at Cerro Pachon (Claver et al. 2004). Because the weak- 
lensing shape measurements will be carried out in r and 
i band images, it is critical for the r and i band imag- 
ing to provide the excellent seeing while still covering the 
required sky area of 20,000 square degrees. From Cerro 
Pachon, the required 20,000 square degrees of survey area 
are viewable above a zenith angle of 45 degrees. At this 
45° limit, the corresponding dispersion in the r and i 
bands is less than ~ 0.3", which elongates the PSF in 
the zenith angle direction by ~ 9% and gives a coherent 
PSF ellipticity of < 0.04 for ~ 0.7" seeing. This amount 
of ellipticity by atmospheric dispersion is small and com- 
parable to the values induced by the optical aberrations 
(after atmospheric turbulence). Therefore, we conclude 
that the LSST design does not require atmospheric dis- 
persion corrector (ADC). 

3. IMPLEMENTATION OF ATMOSPHERIC TURBULENCE 

After the optics and camera are corrected by the wave- 
front sensing, the largest residuals are due to the at- 
mosphere and the focal plane non-flatness. The simu- 
lation of the atmospherically distorted wavefront is the 
first step for the generation of the realistic PSF of LSST. 
A typical short exposure image of a star through atmo- 
spheric turbulence can be viewed as a superposition of a 
number of speckles. The spatial extent of each speckle is 
roughly determined by the diffraction limit of the instru- 
ment whereas the relative centroid shift between these 
speckles reflects the severity of the turbulence mainly 
depending on the amplitude of the low frequencies. To 
mimic these effects numerically, we follow the well-known 
classical approach which regards the atmosphere as con- 
sisting of multiple moving layers each with independent 
phase screen and velocity. Although a single phase screen 
has been shown to provide many characteristics of the 
real environments qualitatively, the wind speed must be 
unrealistically high ~ 100 km s _1 in order to correctly 
represent the time scale of the observed speckle changes. 
In principle, the turbulence affects both the magnitude 
and phase of the wavefront. However, because the phase 
distortion has a much larger effect for non-saturated scin- 
tillation (Fried 1994), in this paper we only consider the 
phase effect on the PSF behavior. 

The power spectrum of the phase screen is usually 
described by the Kolmogorov/von Karman model (Kol- 
mogorov 1941; von Karman 1930; Sasiela 1994): 
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where v, tq, and Lq are the two-dimensional frequency, 
the Fried parameter, and the outer scale, respectively. 
The outer scale Lq is infinite for Kolmogorov screens. 
However, for a large telescope this approximation be- 
comes invalid, leading to an over-prediction of the PSF 
size (e.g., Tokovinin & Travouillon 2006). We set Lq to 
25 m, which is close to the average value at Cerro Pachon 
determined by Ziad et al. (2000) from the measurement 
of the angle- of- arrival fluctuations. 



The realization of the Kolmogorov/von Karman phase 
screen in the spatial domain is obtained by taking the 
inverse Fast Fourier Transform (FFT) of if>(v), whose 
modulus is given by L^/W{y)] L refers to the physical 
length of the screen and the phase of ip{v) is random. For 
the two-dimensional array, where i and j are the indexes 
in the spatial domain, and k and I are the indexes in the 
frequency domain, tp(i,j) (the phase value at i and j) is 
given by (e.g., Carbillet & Riccardi 2010) 
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tion 2 FFTjj/j means that either real or imag- 
art of t 



In equati' 

inary part of the inverse FFT result is used (thus the 
additional factor y2 is present). <f>(k,l) is the uniformly 
distributed random number in the range — n < <fi < ir. 
Because we used FFT, one obvious drawback of the phase 
screen generated in this way is the lack of frequencies 
lower than the value set by the physical dimension of 
the array (i.e., l/L). Lane et al. (1992) suggest a com- 
pensation of this artifact by adding subharmonics to the 
result. However, this problem is not relevant in the cur- 
rent study because the simulation of LSST requires a 
large (> 1 km) phase screen to account for the large field 
of view. 

For the realization of the specific environment for 
LSST, we adopt the results of Ellerbroek (2002), who 
provides the parameters of the atmospheric turbulence 
profile based on the measurements at Cerro Pachon, 
Chile, the site of the current Gemini-South telescope and 
also the future LSST observatory. The atmospheric pro- 
file consists of six layers at altitudes of 0, 2.58, 5.16, 7.73, 
12.89, and 15.46 km as shown in Figure [4] with relative 
weights of 0.652, 0.172, 0.055, 0.025, 0.074, and 0.022, re- 
spectively (Ellerbroek 2002). The Fried parameter (ro) 
is set to 0.16 m, which is reported to be an approximate 
median condition at the location by Ellerbroek (2002); 
for a large telescope with a diameter D » short- 
exposure seeing is inversely proportional to r . Each 
layer is allowed to move independently at a constant ve- 
locity during the integration. The maximum wind ve- 
locity at the highest altitude is kept under ~ 20 m s _1 , 
which determines the minimum time step 0.005 s (corre- 
sponding to ~ 0.5 pixel shift) to satisfy Nyquist sampling 
rate. We choose the dimension of our resulting phase 
screen array to be sufficiently large (6 x 8192 x 8192 or 
6 x 1.3km x 1.3km) in order to keep the moving phase 
screens covering the telescope field of view during the 15 
s integration time. 

The weighted average ip(x,y) of the six phase screens 
at a given moment can be converted to a snapshot of the 
telescope PSF in the absence of the optical aberration 
via the following equations: 



p(x,y) = A{x, y y 2 ^ x >vV x 



PSF = \FFT(p)\ 2 , 
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In equation |3j A(x, y) is the telescope pupil function (a 
mask showing the obscuration) and the phase difference 
i/j(x, y) has the dimension of length. The pixel scale of 
the PSF array given by equation His simply -FA/2, where 
F is the f-ratio. Figure [5] displays the examples of the 
LSST PSFs generated in this way at different integra- 
tion times. The PSF at t — shows the typical instan- 
taneous speckle. In this figure wc do not include either 
the optical aberration or the charge diffusion by CCDs, 
the effects of which we however later add to generate the 
simulated LSST images. As the exposure time increases, 
more speckles are stacked together, which makes the re- 
sulting PSF rounder and the irregular features present 
in the individual speckles more smeared (de Vries et al. 
2007). 

A quantitative study on the impact of the atmospheric 
turbulence on the ellipticity and its spatial correlation 
is needed to support the validity of our simulation here- 
after. As is discussed in £|4j we model LSST PSF varia- 
tion CCD-by-CCD with polynomials. If the anisotropic 
power from the atmosphere within 15 sec exposure turns 
out to be too strong and changes on a very small scale, 
any conventional interpolation scheme with the finite 
number of stars will over-smooth the inherent PSF vari- 
ation. 

We define the ellipticity of PSF as (a—b)/(a+b), where 
a and b are the semi-major and -minor axes of the PSFs, 
respectively. In general, the PSF isophotes change with 
radius. Thus, the measurement somewhat depends on 
the weighting scheme. We measure the ellipticity of the 
PSF using the following quadrupole moments: 

Wij jd 2 ew{e)i(9) ' { ' 

where 1(9) is the pixel intensity at 8, is the center 
of the star, and W(0) is the optimal weight function 
required to suppress the noise in the outskirts. For the 
diffraction-limited PSF of LSST (Figure [8b, we choose a 
Gaussian with a FWHM of ^0.3"for WJ&) whereas for 
the turbulence-limited PSF the FWHM of the Gaussian 
weight function used is ~ 0.7". The quadrupole moments 
are converted to the ellipticity e via: 

e = (e+,e x ) 

_ / Qn — Qii 

VQn + Q22 + 2(QuQ 22 - Qf 2 )V2 ' 
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The magnitude of the stick is proportional to |e| = 

+ £ x 1 ano - the orientation angle is given by 

0.5tan _1 (e x /e + ). 

Figure [6^1 displays the PSF ellipticity variation for the 
15 s exposure within a 4kx4k CCD of LSST by atmo- 
sphere. Because no optical aberration is introduced yet, 
this 4000 "whiskers" show the purely atmospheric contri- 
bution. The size of the sticks represents the magnitude of 
the ellipticity whereas the the orientation of the stick is 
aligned with the position angle of the elongation. A clear 
spatial correlation is visible on a scale of ~1' (~300 pix- 
els). However, the average magnitude is less than 1%, 



and thus the resulting anisotropic power is not strong. 
Figure [6]d shows the ellipticity correlation based on the 
following equation: 

£fc(r) = (et(r') e t (r' + r) ± e x (/) e x (/ + r)), (7) 

where the product e(r') e(r' + r) is performed for each 
pair separated by r, and the e t and e x are the tangen- 
tial and the 45° components, respectively, with respect 
to the line connecting the pair. The amplitude of the 
correlation function decreases with exposure time. From 
the comparison with the contribution due to the tele- 
scope and camera aberration (see Q, we conclude that 
the small scale anisotropic power for the delivered PSF 
is not dominated by the atmospheric turbulence. These 
results have been validated with a series of 15 sec ex- 
posures of the globular cluster NGC2419 on the Subaru 
telescope. 

4. LSST FOCAL PLANE AND IMPACTS ON PSF 

The 64 cm diameter focal plane of LSST will be tiled 
with 189 4kx4k CCDs. As-built heights might vary up 
to ~ 10 microns (peak-to-peak) in a complicated way rel- 
ative to the nominal flat surface. Although this flatness 
deviation is small compared with other cameras in exist- 
ing facilities, the small f-ratio of the instrument makes 
this focal plane flatness variation play a critical role in 
aberration-induced PSF behavior (depth of focus oc f- 
ratio). In Table 1 we summarize the current specification 
of the focal plane error budget originally issued to ven- 
dors, as well as the values used for the current simulation. 
The CCD flatness technology and testing is reviewed in 
Takacs et al. (2006) and Radeka et al. (2009). The 
focal plane errors that we assume for the simulations in 
this paper are larger than the current expectation so that 
the results presented here are conservative. Several ven- 
dors are now exceeding these flatness specifications. A 
realization of the LSST focal plane based on these er- 
ror distributions is displayed in Figure [7J The resulting 
behavior of the PSF ellipticity is illustrated in Figure [8] 
These PSFs are obtained by evaluating optical path dif- 
ference functions across the focal plane in the absence of 
atmospheric turbulence using the ZEMAX software; see 
also Jarvis, Schechter, & Jain (2008) for the discussion on 
the impact of the telescope focus on the PSF ellipticity 

Several features are noteworthy in Figure [HJ First, the 
aberration-induced ellipticity is large. Most of the large 
sticks in the plot exceed ~ 10% ellipticity, approaching 
nearly 30% at the field edges. Because the maximum de- 
viation in this realization is small (~ 10 /im), these large 
values of ellipticity remind us that the focal error toler- 
ance of LSST is indeed narrow; however, we note that 
in the central region (< 1°) the sensitivity of PSF elon- 
gation to height error is somewhat mitigated. Second, 
sharp discontinuities of PSF ellipticity are present where 
the CCD heights also vary abruptly. Third, although it 
may be a bit difficult to recognize this feature in Fig- 
ure [8j a small-scale variation is observed even within a 
CCD mainly due to the tilt and potato chip effects. 

The presence of these smooth, small-scale variation 
and discontinuous changes across CCD gaps are the most 
challenging aspects of the PSF description and modeling 
for LSST. Any attempt to use a single set of polyno- 
mials to characterize the PSF behavior across the en- 
tire focal plane fails because the small-scale variation re- 
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quires impractical, high-order terms in the polynomials, 
and in addition no interpolation scheme can satisfacto- 
rily reproduce the sharp discontinuities across the CCD 
borders. This is the reason that in the current paper 
we perform interpolation CCD-by-CCD to model LSST 
PSFs. Considering the 1-2 cycles of the potato chip ef- 
fect within a CCD, we estimate that 3-4 order polyno- 
mials suffice. The remaining question is what features 
of PSF are interpolated. A simple description of the 
two components of the PSF cllipticity has been used in 
the early weak-lensing analysis (e.g., Valdes, Jarvis, & 
Tyson 1983; Kaiser, Squires, & Broadhurst 1995). How- 
ever, certainly ellipticity alone does not fully character- 
ize a PSF. One obvious missing piece of information is 
the size of the PSF. One can imagine that a 10% el- 
lipticity PSF with FWHM=1" more strongly affects the 
shape of galaxies than a PSF with the same ellipticity 
but with FWHM=0.5". Many authors (e.g., Hoekstra et 
al 1998, Kaiser 2000; Rhodes, Refreiger, & Groth 2001) 
suggested modifications to this earlier method (Kaiser, 
Squires, & Broadhurst 1995) for the treatment of the PSF 
size. Nevertheless, this increased sophistication (i.e., el- 
lipticity with size) is not still sufficient to obtain the ac- 
curacy that the LSST weak-lensing science requires. 

On the other extreme, one might consider interpolating 
pixel values of observed PSFs. The result of this type of 
PSF interpolation is the actual image of the PSF, which 
in principle contains the most complete information at 
any point in the focal plane. However, this scheme faces 
severe practical challenges in the implementation because 
1) except for the pixels near the core, individual PSF 
images are noisy, and 2) the number of needed polyno- 
mials become prohibitively large as one enlarges the di- 
mension of the PSF. Therefore, one might imagine that 
a new PSF description scheme should provide sufficient 
details with both an efficient noise filtering method and 
a compact set of polynomials. A natural choice is the de- 
scription of the PSFs using some basis functions, prefer- 
ably an orthogonal set to avoid degeneracy and thus im- 
prove compactness. Bernstein & Jarvis (2002) and Rc- 
fregier (2003) suggest the use of eigenfunctions of two- 
dimensional quantum harmonic oscillator (also known as 
shapelets) for this purpose. Although the shapelet ap- 
proach is shown to work in many applications, Jee et al. 
(2007) demonstrate that the basis functions are still sub- 
optimal, truncating the PSF profile much earlier than 
desired. Consequently, Jee et al. (2007) proposed to 
derive the basis functions from the data themselves via 
principal component analysis (PCA). 

Lupton et al. (2001) first mentioned the use of PCA for 
the modeling of the Sloan Digital Sky Survey (SDSS). Al- 
though the details of the implementation were not given, 
the method is similar to that of Jee et al. (2007) in 
that PCA is employed to derive the optimal basis func- 
tions from the data themselves. On the other hand, the 
use of PCA by Jarvis & Jain (2004) and Schrabback et 
al. (2010) for the description of the PSF variation is 
somewhat different from the method of Jee et al. (2007) 
because they use PCA to characterize the pattern of the 
variation (e.g., coefficients or cllipticity components), not 
to capture the pixelated PSF features. 

Our PCA approach has been so far applied to weak- 
lensing studies only with Hubble Space Telescope images 
(e.g., Jee et al. 2009). In the current investigation, we 



adopt this PCA scheme to characterize/model the LSST 
PSFs. As our PCA algorithm is independent of the spec- 
ification of the instrument, the implementation in Jee et 
al. (2007) can be used with little modifications. The 
technical details of the PCA method will be reviewed in 
the following section. 

5. PRINCIPAL COMPONENT ANALYSIS 

PCA is the mathematical procedure for analyzing mul- 
tivariate data by finding and utilizing a new orthogonal 
coordinate system, which substantially reduces the num- 
ber of axes relative to the original number of variables 
for the characterization of the data. Inevitably, the pro- 
cedure involves information loss. However, in the typical 
data set where the noise is significant, this reduction of 
the information in fact works as filtering, enabling a com- 
pact description of the data without significantly harm- 
ing its integrity. In the current paper, we explain the 
procedure using our specific data set, i.e., series of noisy 
two-dimensional PSF images. A similar account can be 
also found in Jee et al. (2007). 

Let us assume that we have N observed stars in a 
4kx4k CCD image with random positions and bright- 
ness. The two-dimensional array of each PSF image can 
be rearranged to form a one-dimensional vector. How 
we exactly scramble the pixels of the PSF image is not 
important in the subsequent analysis as long as we have 
the ability to re-order the pixels to recover the original 
two-dimensional image. However, it is critical that each 
PSF image is carefully background-subtracted and nor- 
malized prior to the rearrangement. Then, the N one- 
dimensional vectors are stacked to comprise a matrix. If 
each PSF image contains M pixels, the result isaiVxM 
matrix. The next step is to evaluate mean values along 
the columns, and subtract this lxM mean vector from 
each row of the matrix. In practice, each row of the 
N x M matrix contains star images contaminated by 
neighboring objects, cosmic rays, bad pixels, etc. There- 
fore, the mean values are preferably evaluated iteratively 
with some outlier rejection. Also, the rejected pixel val- 
ues in the N x M matrix should be replaced with the 
mean values. Because we need this mean vector to re- 
construct the PSF later, this mean vector (or mean PSF 
rearranged in one-dimension) must be saved. 

Now, we are ready to find a set of P (< MIN{N, M}) 
orthogonal vectors that compactly describe the N x M 
matrix. One way to derive such orthogonal vectors is 
Singular Value Decomposition (SVD; Press et al. 1992). 
If we call the N x M matrix S, the matrix S can be 
decomposed in the following way: 

S = UWV T , (8) 

where UisaniVxM column-orthogonal matrix, W is an 
M x M diagonal matrix, and V is an M x M orthogonal 
matrix. The diagonal elements of W are called singular 
values. 

Therefore, the matrix component Sjj can be expressed 

as 

M 

Si 3 = ^2w k U lk V jk -, (9) 
fe=i 

which implies that, if some singular values w k are rela- 
tively small, we can substitute zeros for those w k 's with- 



G 



Jee et al. 



out greatly compromising the numerical accuracy of the 
matrix evaluation. In other words, the reconstruction of 
S is possible with fewer columns in both U and V, and 
thus the procedure provides a substantial compression of 
the information in the original matrix S, where photon 
and detector noise is a significant part. The remaining 
columns (principal components) of V form an orthogo- 
nal basis, and in this analysis we refer to these eigen- 
functions as eigenPSFs. We display examples of these 
eigenPSFs in Figure [9j These eigenPSFs in this fi gurc 
are obtained from the PCA of 4000 simulated stars on 
the entire focal plane in the absence of telescope aberra- 
tion and CCD height fluctuation, but in the presence of 
the atmospheric turbulence. The eigenPSFs are ranked 
according to the size of their eigenvalues with n = rep- 
resenting the largest. Note the remarkable resemblance 
of the first several eigenPSFs to the basis functions in 
the shapelet approach. However, also note the asymme- 
try of the eigenPSFs due to the characteristic pattern 
of the specific data set. For increasing n, more high- 
frequency features are evident. Above a certain limit 
n = n c , the resulting eigenPSFs start to reflect the noise 
features without contributing to the PSF reconstruction. 
Consequently, we need only a truncated set of eigenPSFs 
to describe the behavior of the PSF. 

This minimal number of eigenPSFs can be determined 
by examining the eigenvalues of each eigenPSF. In gen- 
eral, there is a conspicuous "kink" at a certain n = n c , 
beyond which the decrease of the eigenvalue suddenly 
slows down (e.g., Figure 5 of Jee et al. 2007). We find 
that ~ 20 eigenPSFs are sufficient to describe the details 
for both space- and ground-based telescope PSFs. 

After the eigenPSFs are obtained, the k th star used for 
the derivation is decomposed as 

Ckihj) = £ a kn P n (i,j)+T(i,j), (10) 

where Ck(i,j) is the normalized pixel value of the k th 
star at the pixel coordinate (i, j), P n is the n th eigenPSF, 
cifc„ is the projection of the k th star in P n , and T is the 
mean PSF. Because P n 's are orthogonal to one another, 
one does not need to determine a,k n through some % 2 
minimization here. 

Now the PSF at any arbitrary location can be evalu- 
ated through polynomial interpolation of cifc n 's: 

l+m<N 

a„ = ^2 d ni m x l y m (11) 

Z,m=0 

where x and y are the coordinate of the position, at which 
one wants to know the PSF, d n i m is the coefficient of the 
term x l y m for the n th eigenPSF, and N is the maximum 
order of the polynomial. It is important to determine 
d n im without the application of the S/N weight of the 
star because we are interested in the spatial variation of 
the PSF. 

6. SIMULATION OF LSST IMAGES AND PSF 
RECONSTRUCTION 

For the robust simulation of the weak-lensing shear 
extraction using LSST, it is crucial that the simulated 
image contains astronomical objects with realistic photo- 
metric, morphological, and statistical properties, as well 



as the instrumental artifacts resulting from the known 
characteristics of the telescope, camera, and the environ- 
ment. Specifically, the requirements include our ability 
to correctly represent: 

1. the number density and luminosity function of un- 
saturated high S/N Galactic stars, 

2. the number density, luminosity function, size, and 
morphology of galaxies, 

3. position-dependent PSF variations, 

4. geometric distortion, and 

5. atmospheric turbulence and dispersion. 

The first requirement is needed because the number of 
available high S /N stars determines how well we can mea- 
sure the behavior of the PSFs. In particular, we need to 
investigate if the nominal 15 s exposure image provides 
sufficient number of usable stars within a CCD so as to 
correctly sample the position-dependent PSF variation. 

The next point relates to the issue of both statisti- 
cal and systematic errors of the shear measurement from 
galaxy ellipticity. The number of usable galaxies for test- 
ing lensing analysis affects the statistical accuracy of the 
estimate of shear and our determination of mass. In ad- 
dition, one important source of bias in many ellipticity 
measurement techniques is the under-fitting of galaxy 
shapes (i.e., the model galaxy profile assumed is differ- 
ent from the real profile). Since many z > 1 galaxies are 
dominated by so-called faint blue galaxy (FBG) popula- 
tion incapable of being described by analytic functions, 
it is important that the image contain a realistic fraction 
of these irregular galaxies. The extremely high sensi- 
tivity of the LSST PSF to the focal plane errors makes 
the position-dependent PSF one of the most critical fea- 
tures to be simulated. We need to generate the realistic 
PSF and convolve the galaxy image with this position- 
dependent PSF in a time-efficient way. 

Although expected to be very small, the geometric dis- 
tortion by the telescope and atmosphere can correlate el- 
lipticity of galaxies and mimic gravitational lensing sig- 
nal similar to the way that uncorrected PSF can cause 
systematic alignment of galaxy ellipticity. It is essential 
to estimate the level of accuracy in rectifying the raw 
images and the effects on shear measurement. 

The atmospheric turbulence and dispersion are dom- 
inant sources of time-dependent PSF variation. Un- 
like the telescope and camera optical aberration, the at- 
mospheric turbulence is believed to circularize the PSF 
rather than to introduce significant anisotropy whereas 
the atmospheric dispersion obviously elongates the PSF 
in the direction of the zenith angle. It goes without say- 
ing that these atmospheric effects should be properly in- 
cluded in the simulation. 

Our simulation of the LSST shear test is an on-going 
project, and the final goal is to carry out an end-to-end 
simulation from the realization of the LSST image from 
an input cosmology to the full recovery of the galaxy 
shear signal and the cosmological parameters. Hence, it 
is our ultimate aim to rely on the final simulator not 
only to implement the above requirements but also to 
address other non-shear measurement issues such as dif- 
fuse components, cosmic rays, hot/warm/bad pixels, etc. 
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This end-to-end simulation is being done by The LSST 
Simulations Team (Connolly et al. 2010). 

Here, as a pilot study, we focus on generating realistic 
LSST PSFs and recovering these input PSFs from sparse 
distribution of the stars in the presence of background 
galaxies. Therefore, some simplifications are made where 
they do not significantly compromise the integrity of our 
simulation relevant to the focal plane flatness issue. 

6.1. Generation of Input PSF Model 

First, we generate the coordinates of the model stars 
that are densely distributed across the focal plane. Be- 
cause we will evaluate the fiducial PSFs at these loca- 
tions and interpolate the results to describe the PSF on 
every pixel, the density of these stars should be high 
enough to accommodate the expected small scale vari- 
ation. We assigned about 800 stars to each CCD and 
verify that these 800 stars sample with ample resolu- 
tion the spatial variation, which is dominantly modu- 
lated by potato-chip effects (1-2 cycles per CCD). Then, 
we evaluate the diffraction-limited PSF at the location 
of those stars using ZEMAX. Strictly speaking, in this 
step we save phase data rather than PSFs because the 
optical phase data should be later combined with the at- 
mospheric phase screen to generate the final PSF. We 
use the simulated focal plane height fluctuation shown 
in Figure [7] in order to provide the focus value for each 
star. In the current investigation, our main interest is the 
ability to model/recover the PSF in the r band. Thus, 
we simulate the broadband effect by obtaining the data 
at five uniformly spaced wavelengths between 5450-7030 
A. Next, the atmospheric phase screen is evaluated at 
the location of each star. This atmospheric phase infor- 
mation is merged with the optical data to construct the 
next stage PSF. Finally, the result is further convolved 
twice to mimic the effects of the charge diffusion (<~ 0.4 
pixel rms) within a CCD and the atmospheric dispersion 
(close to zenith). 

We derive sets of eigenPSFs separately for each 
CCD. Alternatively, one can consider finding one set 
of eigenPSFs utilizing all PSFs on the focal plane. Al- 
though this latter scheme may reduce statistical errors, 
it appears that the resulting eigenPSFs are subopti- 
mal, somewhat under-predicting the systematic pattern 
unique in each CCD. We choose to keep the most sig- 
nifincant 20 eigenPSFs, which accounts for more than 
99% of the total variance. Now each PSF is described by 
20 coefficients each representing the amplitude in each 
eigenPSF direction. We fit 4th order polynomials to the 
coefficients of all stars to describe the PSF at any arbi- 
trary location. 

6.2. LSST Galaxy Image Generation 

One ideal way to generate galaxy images for a given 
cosmology is to make a catalog using large numerical 
simulation results (e.g., Millennium Simulation, Springel 
et al. 2005) with ray-tracing technique, to assign proper 
spectral energy distribution (SED) and morphology to 
individual galaxies, and to convert the catalog to as- 
tronomical images. However, this full scale end-to-end 
simulation (Connolly et al. 2010) is beyond the scope 
of the current project, and here we instead utilize Hub- 
ble Space Telescope (HST) Ultra Deep Field (UDF) im- 
ages to sample galaxies. The HST UDF (Beckwith et 



al. 2003) is a 412-orbit HST Cycle 12 program to im- 
age a single field with the Wide Field Camera (WFC) of 
the Advanced Camera for Surveys (ACS) in four filters: 
F435W, F606W, F775W, and F850LP. The total expo- 
sure time of the F775W image is ~ 350 ks, and the 10 
a limiting magnitude is 29 ABmag. The depth of the 
UDF image is sufficient for the generation of the galax- 
ies that will be imaged by LSST, which reaches r ~ 27.5 
mag (~ 5c) in the final depth. Hence, virtually, the UDF 
data provide noiseless, seeing-free galaxies (the size of the 
ACS PSF is less than that of one LSST pixel) for a 15 s 
exposure. As the field of view of the ACS is ~ 3'x3', even 
smaller than the area covered by one LSST ~ 14' x 14' 
CCD, it is inevitable that the same UDF galaxies must 
appear multiple times across the focal plane. However, 
this is not likely to compromise the accuracy of the cur- 
rent weak-lensing simulation, where we do not attempt 
to recover any input cosmology. 

The high-level science product of the UDF is publicly 
availably and we use the F775W image to sample galax- 
ies. Objects were detected via SExtractor (Bertin & 
Arnouts 1996) by looking for 8 continuous pixels above 
the sky rms. Postage stamp images were created by crop- 
ping a square region centered on the object with the side 
five times the semimajor axis of the SExtractor estimate. 
We used the segmentation map to check if any pixels be- 
longing to other objects were present in the square. If 
so, we filled these pixels with the randomized background 
whose statistics were tuned to match the SExtractor es- 
timate. Any objects whose half-light radii were less than 
the upper limit of the HST PSF value were discarded. 

The next procedure is to plant these postage stamp 
UDF galaxy images onto each 4k x 4k CCD image. We 
looped over the list of the postage stamp images sequen- 
tially but randomizing the location and the orientation. 
The loop ended when the desired number density was 
reached. A large fraction of the pixels were the result of 
the superposition of many different postage stamp im- 
ages. Thus, it is apparent that the noise level of these 
pixels vary slightly depending on the number of overlaps. 
However, because we later introduce much larger addi- 
tional noise to simulate the 15 s exposure data, this noise 
variation is negligible. 

A stellar count study based on the SDSS catalog (Juric 
et al. 2008) shows that a single 4k x 4k CCD image 
(~ 14" x 14") in 15 s exposure at the Galactic north 
pole is expected to contain at least ~ 140 high S /N stars 
(> 30<r), which can be used to model the PSF. However, 
in the current simulation we conservatively distributed 
40 point sources per CCD in such a way that the S /N of 
those stars remain above 30 in 15 s exposure, assuming 
the following differential star count: 4>(m r ) cx lO - 27 ™ 1 -; 
we find that PSF reconstruction with different magni- 
tude distribution does not impact the quality of the PSF 
reconstruction as long as all the stars used have S/N> 30. 

Prior to the noise generation, the "perfect" 4k x 4k 
image must be convolved with the position-dependent 
PSF of the CCD. Because our PSF is represented by 20 
coefficients and 20 basis functions, it provides us with a 
convenient way to implement the spatially varying con- 
volution. To begin with, we convolve a 4kx4k seeing-free 
astronomical image with 20 eigenPSFs via FFT. The re- 

3 http:/ /archive. stsci.edu/prepds/udf/ 
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suits of this operation are 20 4kx4k images. Then, using 
the resulting 4th order polynomial fitting we generate 
another 20 4kx4k images with each pixel containing the 
coefficient (amplitude) of the eigenPSF. Next, we mul- 
tiply each eigenPSF-convolved image with the match- 
ing coefficient image, obtaining 20 new 4kx 4k images. 
These 20 images are simply stacked together to produce 
one 4kx4k image. This image is not yet the final im- 
age that we desire to generate because the 20 eigenPSFs 
are derived from the mean-subtracted PSFs. Therefore, 
we need to convolve the seeing-free 4kx4k image with 
the mean PSF, and add the result to the former 4kx4k 
image to yield the final 4kx4k LSST image. Figure [TO] 
displays a 100" x 100" cutout from one of the 189 4k x 4k 
images. The noise-free (except for the existing low-level 
noise in the ACS UDF image) LSST-seeing-convolved im- 
age is shown in the left panel. This "ultra deep" image is 
about 2 mag deeper than the median depth LSST image 
(middle) after the completion of the nominal 10 year mis- 
sion with 200 exposures per sky patch. The nominal 15 
s LSST single exposure simulation is shown in the right 
panel. 

7. SIMULATED IMAGE ANALYSIS AND RESULTS 
7.1. PSF Reconstruction 

The procedure to apply PCA to the simulated LSST 
image containing galaxies and stars is similar to the one 
discussed in the case of the input PSF modeling. The 
difference is that now we have far fewer stars, realistic 
noise, and neighboring objects. As a matter of course, 
this will provide the reconstructed PSF with fewer de- 
grees of freedom in the description of the spatial varia- 
tion. Hence, if the intrinsic PSF variation within a CCD 
is more complicated than what the reconstructed model 
can predict, the weak-lensing science with LSST will be 
substantially compromised. In this section, we present 
the detailed comparison between the reconstructed PSF 
and the input PSF in order to address this critical issue. 

Astronomical objects in the simulated 189 4kx4k im- 
ages were detected with SExtractor. Stars were selected 
by utilizing half light radius versus magnitude relation as 
performed in typical lensing analysis. The stellar locus 
is well defined for bright stars. In addition, using ellip- 
ticity, comparison of half-light radius with FWHM mea- 
surement, and source extraction history (e.g., presence 
of close neighbor, edge clipping, etc.), the contamination 
becomes negligible in the current 10 sq. degree simula- 
tion. The cutout size is chosen to be 31pixel x 31pixel 
(6.2" x 6.2") after the star is sub-pixel shifted to place 
the centroid on a center of pixel. The ellipticity of the 
selected stars are displayed in Figure 11 We overlay the 
"whiskers" over the focal plane image that we use for 
the input PSF generation. The pattern of the PSF vari- 
ation resembles the case displayed in Figure [8] where no 
atmospheric effect is included. However, the reduction 
of the magnitude is apparent due to the circularization 
effect of the atmospheric turbulence. We observe that no 
significant anisotropy is introduced by the moving phase 
screens. The PSF elongation in the zenith angle direction 
is very small here because we choose a telescope pointing 
to be close to the zenith. 

We derive eigenPSFs for each CCD separately and use 
3rd order polynomials to describe the spatial variation. 



For the visualization of how well the reconstructed model 
represents the observed pattern, we display the resi dual 



ellipticity distribution across the focal plane in Figure 12 



The residual in Figure [12] is the difference in ellipticity 
(e + ,e x ) of the observed (simulated) stars and the model 
PSF evaluated at the location of the stars. Consequently, 
any correct model sh ould not only reduce the size of the 
"whiskers" in Figure 11 but also make the correlation 



of their size and the orienta tion decrease sig nifi cantly. 
Visual comparison of Figure 1X2] with Figure [TT] shows 
qualitatively that our interpolated PSF at thelocation 
of the stars closely mimics the pattern of the observed 
PSFs. One method to quantitatively examine the reduc- 
tion of the whiskers is to compare the ellipticit y c om- 
ponent distribution in the two cases. In Figure 13 the 
blue dots represent the ellipticity components of the ob- 
served stars whereas the red dots indicate those of the 
residuals. Note that the residual ellipticity distribution 
is much more compact and is precisely centered on the 
origin (e^,e x ) = (0,0). 

Figure [14] displays the ellipticity correlation. The spa- 
tial correlation of the observed ellipticity (before PCA 
correction) has an amplitude of 3 x 10~ 5 . For a sanity 
check, we compared this amplitude with the values that 
we obtain from real observations, and found that they 
agree well. For example, the amplitude of the intrin- 
sic PSF correlation before correction is ~ 4 x 10~ 5 and 
~ 3x 10 -5 in the Deep Lens Survey (DLS; Wittman et al. 
2002) field from the 4m Kitt Peak telescope and 8m Sub- 
aru telescope images, respectively; one should use cau- 
tion when comparing these values to other published re- 
sults because of different ellipticity definitions and weight 
functions. 

Also displayed in Figure [Ti] is the ellipticity correlation 
function of the purely atmospheric PSF (Figure[6]), which 
is generated by "turning off" optical aberrations. These 
PSFs are free from interpolation errors and thus help us 
to assess the amount of the atmospheric contribution to 
the PSF anisotropy. In particular, the comparison en- 
ables us to verify that ~ 800 stars per CCD is sufficient 
to model the PSF pattern within a CCD. If the atmo- 
spheric PSF possessed a significant power on a small scale 
r < 1', those ~ 800 stars could not adequately sample 
the PSF variation. 

Obviously, the amplitude 3x 10 -5 of the ellipticity cor- 
relation for the observed PSF (blue) is higher than the 
expected cosmic shear signal for the A CDM universe at 
large angles. This implies that without PSF correction 
the galaxy ellipticity correlation will be dominated by 
the PSF systematics rather than by the intrinsic lens- 
ing signal. It is reassuring that the correlation function 
of the residual ellipticity components has an amplitude 
of ~ 10 -8 , which is more than 3 orders of magnitude 
reduction in correlation and also more than 2 orders of 
magnitude below the anticipated cosmic shear lensing 
signal (see Rowe 2010 for the discussion and use of the 
residual ellipticity correlation in detail). We find similar 
improvement when we apply the PCA PSF corrections 
to the DLS data. 

It is important to remember that the above PSF cor- 
relation is measured on a single 15 s exposure image. In 
fact, LSST will visit the same patch of the sky roughly 
200 times in each filter used for weak-lensing with two 15 
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s exposures. Before each exposure the telescope optics is 
updated from the wavefront sensing system. The pupil 
rotation, sky rotation, and x — y dithers will mix different 
systematic PSF patterns and reduce the intrinsic ellip- 
ticity of the PSFs in the final stack. In addition, this will 
also remove the spatial correlation that we saw within a 
single CCD even before we apply PSF corrections. 

Although one can come up with an optimized opera- 
tion scheme, which takes full advantage of hundreds of 
visits, we observe that even a simple random field rota- 
tion can dramatically reduce the intrinsic PSF ellipticity 
correlation. Figure [15] schematically displays this sim- 
ple observing pattern and the resulting PSF ellipticity 
distribution. We illustrate an over-simplified observing 
scenario, where a same patch of the sky is visited 100 
times with randomized roll angles. 

In order to create the stack image out of these 100 
frames, we arbitrarily choose one frame as a reference and 
rotate the rest of the frames for the alignment. Stacking 
rotated images are highly non-trivial in the presence of 
anisotropy in pixel boundaries, atmospheric dispersion, 
satellite trails, cosmic rays, detector response variation, 
internally scattered light/ghosts, etcrlNumerous schemes 
exist for fractional pixel handing. VVe use bi-cubic in- 
terpolation (Press et al. 1992), which closely mimics 
the ideal sine interpolation or its variations (e.g., Lanc- 
zos kernels), for both stacking images and constructing 
PSFs; the merits of using Lanczos kernels are discussed 
in the literature (e.g., Jee et al. 2007; Bertin 2010). The 
atmospheric dispersion is very small in the current sim- 
ulation because we simulated pointings near the zenith. 
Thus, the dispersion correction is omitted in the current 
simulation. However, we think this atmospheric disper- 
sion and other aforementioned real world problems must 
be handled with care because any small bias may appear 
as second-order systematics after analysing many billion 
galaxies in the 40,000 sq. degree area. These issues will 
be discussed in our future publications. 

We use the single realization of the focal plane (shown 
in Figure [7| in the above multi-epoch simulation al- 
though we randomize the phase of the atmospheric tur- 
bulence screen for each epoch. It is important to realize 
that this simulation setup for the telescope and camera 
closely reflects the planned operation scheme of LSST. 
The LSST camera will maintain its fixed focal plane con- 
figuration throughout the mission, and the mirror de- 
formation and alignment changes will be corrected with 
actuators to < 0.2/xm accuracy (Manuel et al. 2010). 
That is, the difference in aberration for different visits 
is small compared to the aberration introduced by the 
focal plane non-flatness. In principle, this small wave- 
front reconstruction error may slightly alter the result- 
ing telescope aberration patterns exposure by exposure. 
However, the wavefront reconstruction error in different 
exposures are random, and stacking many exposures will 
reduce the correlation, circularizing the final PSFs to our 
advantage. 



In Figure 15 the observed PSF ellipticity after stacking 
(upper right) is much smaller than what is seen in indi- 

4 While beyond the scope of this paper, our "Multi-Fit" ap- 
proach (Tyson et al. 2008) to galaxy photometry and shape mea- 
surement fits a galaxy morphological model to all the exposures, 
and thus can avoid the propagation of these issues into the LSST 
weak-lensing science. 



vidual single 15 s observations (upper left). The lower left 
panel shows the reconstructed PSF by first performing 
PCA on each exposure and then stacking all model PSFs 
after correct application of field rotation; the rotated 
PSFs are created by re-sampling with bi-cubic interpo- 
lation, which is identical to the procedure in the image 
stacking above. Stacking different PSF patterns reduces 
the correlation in observed PSF (Figure 16 ) However, it 
does not bring down the amplitude below 10 -7 over the 
entire range of angular scales. This is because our toy 
model observing scheme without x — y dithers does not 
mix different PSFs efficiently. 

7.2. Null Test with Galaxies 

The high-fidelity reconstruct ion o f the observed PSF 
through PCA demonstrated in §7.1| hints at the prospect 
that the resulting galaxy shear measurement can be car- 
ried out with similarly small systematic errors. How- 
ever, there remain several issues. First, anisotropic PSFs 
change the pixel noise properties also anisotropically, and 
thus even with the knowledge of the perfect PSF the 
resulting galaxy shape measurement can be still biased 
along the elongation of the PSFs. Related to this point 
is the centroid bias (e.g., Kaiser 2000; Bernstein & Jarvis 
2002), which causes the object centroid to be more un- 
certain in the direction of the PSF elongation. Second, 
our comparison between the model and the observation 
is limited to the PSFs at the location of the stars. If 
the intrinsic spatial variation of the PSF has higher fre- 
quency than the order of the polynomial (3rd order here), 
it is possible that galaxy shape measurements would be 
carried out with suboptimal PSFs. Third, apart from 
PSF modeling, the difference in galaxy radial profile be- 
tween assumption and reality can cause biases in shear 
extraction as suggested by Bernstein (2010). Fourth, it 
is not yet clear how to handle the pixelization effect opti- 
mally. Because photons are collected in finite CCD pix- 
els, any shape measurement algorithm needs to consider 
this procedure to avoid bias. Similar to the convolution 
by seeing, the pixelization also somewhat destroys the 
information, especially at the core of object, where the 
radial profile is most steep. [However, the low surface 
brightness limit of each exposure allows shape measure- 
ment at large radius.] Fifth, object detection from real 
astronomical images involves many subtleties such as de- 
blending, spurious object identification, etc. The current 
algorithms need to be improved substantially to minimize 
the contamination from these undesired detections or at 
least to quantify the fraction more precisely in order to 
estimate the dilution of the lensing signal. 

In terms of shear extraction, the above concerns can 
be largely divided into two issues: the removal of the 
anisotropic bias and the calibration bias. For conve- 
nience, many authors (Heymans et al. 2006; Massey et 
al. 2007; Bridle et al. 2010) formulate the two issues 
with the following expression: 



7 — 7 = rrvy + c, 



(12) 



where m and c are the calibration and anisotropic biases, 
respectively, and 7 and gamma* are the observed and 
true shears, respectively. Despite the difficulties men- 
tioned above, as a last resort the calibration bias or mul- 
tiplicative factor m can be determined through extensive 
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image simulations. We prefer a method which provides 
self-contained calibration factor and relies on the image 
simulation result only as a confirmation tool. Within the 
scope of this paper, we address the second issue, i.e., the 
removal of anisotropic bias, by performing "null" test on 
the simulated LSST images. In particular, we focus on 
how well a simple field rotation can remove this additive 
factor. 

For the measurement of the galaxy ellipticity, we fit 
PSF-convolved elliptical Gaussians (e.g., Jee et al. 2009) 
to galaxies. As mentioned in §7.1| the PSF on the fi- 
nal stack is constructed from the superposition of the 
100 PSFs from all visits. We discard objects whose half- 
light radius is less than 1.2 times the median value of 
the stars or ellipticity measurement errors greater than 
0.2. The resulting number of usable galaxies is ~ 40 per 
sq. arcmin for the depth at the completion of the nom- 
inal LSST survey, which meets the weak-lcnsing science 
requirement (LSST Science Collaborations et al. 2009). 
The exact value of the selection criteria should be deter- 
mined in conjunction with the study of shear calibration, 
and thus will be the subject of our future publications. 

As no shear is introduced, the galaxy correlation func- 
tion should possess amplitudes consistent with the shot 
noise. One should remember that since the intrinsic el- 
lipticity of individual galaxies are at least a factor of 10 
larger than that of PSFs, the correlation function within 
a single LSST field will not approach the level of the 
residual PS F co rrelation function. 

In Figure [TTJ we present this "null" test result with 
galaxy ellipticity for the 10 sq. degree field. If anisotropic 
PSF effects are corrected, the correlation amplitude 
should simply reflect the statistical galaxy shot noise as 
observed. As the number of exposures increases, the cor- 
relation amplitude goes down because the residual PSFs 
become more uncorrelated (systematic noise) and also 
the increased depth gives higher number of galaxies us- 
able for shear measurements (statistical shot noise) . The 
number density of galaxies with sufficiently high S /N for 
ellipticity measurement for a single, 10, and 100 expo- 
sures are 10, 24, and 42 galaxies per sq. arcmin, respec- 
tively. We find that the mean amplitudes of the galaxy 
ellipticity correlation functions are closely proportional 
to the inverse of these number densities (Figure 18 ). This 
tight relation implies that the level of the systematic er- 
ror for the 100 exposure data set is much less than that 
of the statistical error (< 10~ 7 ). The expected signal 
for a ACDM cosmology corresponds to a galaxy elliptic- 
ity correlation > 10~ 6 over this range of angular scale. 
The exact level of the systematic floor can be estimated 
only by extending the area of the simulation by a large 
amount. For the case that the amplitude of the PSF ellip- 
ticity correlation function ~ 10~ 8 in Figure Il6j happens 
to be the fundamental limit, we estimate the extrapo- 
lated correlation function for the 20,000 sq. degree area 
survey (cyan line in Figure 17 1. 



8. CONCLUSIONS AND SUMMARY 

We have investigated how the atmospheric turbulence 
and the telescope aberrations will impact the PSF of the 
active optics LSST in the context of weak lensing sci- 
ence. The atmospheric turbulence simulated by six mov- 
ing Kolmogorov phase screens mainly contributes to the 
broadening of the diffraction-limited PSF (that would 



have been obtained in the absence of the atmosphere) 
without significantly inducing additional anisotropy. The 
FWHM - 0.7" of the simulated PSF nicely agrees with 
the mean value determined by the LSST seeing monitor- 
ing program at Cerro Pachon. The instrumental aberra- 
tions place potential challenges to precision weak-lensing 
analysis with LSST because of the fast f-ratio of the op- 
tics. Our simulation shows that the resulting shallow fo- 
cal depth creates abrupt changes in PSF pattern across 
CCD borders. As the 3.5° x 3.5° focal plane will be tiled 
with 189 CCDs each having its own surface distortion 
(i.e., potato-chip effect), the PSF pattern on the entire 
focal plane is very complex. 

However, we demonstrate that the complexity of the 
PSF pattern over the entire focal plane can be nicely de- 
scribed by interpolation with basis functions derived by 
PCA. In order to account for the discontinuity across the 
CCD borders with polynomials, the basis functions and 
the polynomial coefficients are determined separately for 
each CCD. This PSF recovery test is carried out with 
simulated LSST images containing realistic galaxies, as 
well as stars in the presence of photon noise and neigh- 
boring objects. The residual PSF ellipticity correlation 
function has amplitudes of 10 -8 over the r = 0.2 ~ 3.0° 
range, significantly lower than the amplitudes of the cos- 
mic shear and the raw PSF ellipticity correlation. 

It is critical to remember that in actual LSST opera- 
tions the same patch of sky will be visited multiple times 
with different combinations of pupil rotation, sky rota- 
tion, and x — y dither. Hence, this will substantially 
reduce our burden in modeling the PSF for individual 
exposures because the mix of different PSF pattens will 
circularize the PSF in the stacked image and also de- 
creases the spatial correlation. To illustrate how much 
the stacking procedure can improve the systematics, we 
simulate a toy case, where 100 visits are randomly ro- 
tated around the center of the field. The reduction in 
both ellipticity and correlation prior to any PSF model- 
ing and correction is dramatic even in this simplistic sce- 
nario. The "null" test with galaxies in the stacked image 
indicates that the amplitude of the galaxy shear correla- 
tion is consistent with the level of the shot noise, indi- 
cating that the systematic error is not detectable within 
a single ~ 10 sq. degree field with 100 exposures. 

The PSF anisotropy correction (additive factor in shear 
recovery) is one of the most poorly investigated effects, 
and the literature disagrees on the best approaches. Our 
current simulation shows that the anisotropic PSF ef- 
fects can be sufficiently corrected with the current de- 
sign of LSST. This improvement is enabled by both a 
new PCA algorithm and the large etendue and active 
optics of LSST capable of turning systematics into statis- 
tics. In a forthcoming paper, we will extend the current 
simulation with the addition of input galaxy shear and 
discuss the accuracy of the shear recovery with existing 
algorithms. 
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Fig. 1. — LSST optical design. The telescope design (left) is a modified Paul-Baker three mirror system, which produces a uniform image 
quality across a large field of view. The relative positions of the primary and the tertiary mirrors were adjusted during the design process 
so that their surfaces meet with no axial discontinuity at a cusp, allowing the Primary and the Tertiary to be fabricated from a single 
substrate. The 3.4-m convex secondary mirror has a 1.8 m inner opening, through which the LSST camera is inserted at the focal plane. 
The focal plane (right) is tiled with 189 4kx4k CCDs. There are four wavefront sensors and eight guide sensors, located in the four corners 
of the CCD array. All optics and the camera are corrected with active optics wavefront curvature sensing. The red solid circle represents 
the 3.5° diameter field of view boundary. 
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FlG. 2. — Encircled energy of LSST PSF. Here we estimate the diffraction-limited encircled energy for i filter. The left panel shows the 
encircled energy as a function of radius for the four locations on the focal plane. On the right panel the 80% encircled energy radius as a 
function of the focal plane location (distance from the principal axis) is plotted. The inset images display the snapshot of the PSFs at 0° 
(center), 0.5°, 1°, 1.2°, and 1.75° (edge). Note the remarkably small variation of the PSF size across the entire focal plane. 
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Fig. 3. — Simulated Atmospheric dispersion at the proposed LSST site. We used the dispersion models summarized in Filippenko (1992) 
assuming / = 8 mm Hg (water vapor pressure), T = 5°C (atmosphere temperature), and P = 520 mm Hg (atmospheric pressure), which 
are the approximate average conditions at Cerro Pachon (Claver et al. 2004). 




Fig. 4.— 



Six layers of Kolmogorov/von Karman phase screens used for the atmospheric turbulence model. 
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Fig. 5. — Example of time evolution of simulated LSST PSF due to atmospheric turbulence at the zenith. The displayed PSFs are for 
the monochromatic beam at A = 6000A. We do not include the optical aberration or the charge diffusion by CCDs in this simulation. 
As the exposure time increases, more speckles are stacked together, which makes the resulting PSF rounder and the surface brightness 
irregularities decrease. The FWHM at t ~ 15s is ~ 0.5". 
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Fig. 6. — Impacts of turbulence on PSF. (a) PSF ellipticity (spatial) variation for the 15 s exposure within a 4kx4k CCD of LSST due to 
atmosphere. Because no optical aberration is introduced yet, this plot shows purely atmospheric contribution. A clear spatial correlation 
is visible on a scale of ~1' (~300 pixels). However, the average magnitude is less than 1%, and thus the resulting anisotropy does not 
contribute significantly to the final PSF. (b) Ellipticity correlation for different exposure times. The amplitude of the correlation function 
decreases with exposure time because the telescope aperture sees higher number of (effectively) uncorrclated phase screens. The solid and 
dashed lines represent and £_ , respectively. 
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Fig. 7. — Simulated focal plane of LSST. We used the CCD assembly /fabrication specification in Table 1 to generate the LSST focal 
plane tiled by 189 CCDs. 
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Fig. 8. — Impact of focal plane CCD height variation on PSF. We use the ZEMAX software to obtain diffraction-limited PSFs in the 
absence of atmospheric turbulence. The encircled stick in the lower-left corner represents 10% ellipticity as illustrated. Note both the 
magnitude of ellipticity induced by the focal plane error (and aberration) and the abrupt changes across the CCD borders. The large circle 
shows the 3.5 m diameter field of view. The sticks outside this circle arc for illustrative purpose only to demonstrate that the aberration 
degrades severely beyond the 1.75 ° boundary. 
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Fig. 9. — PCA cigcnPSFs. The eigenPSFs were ranked by their eigenvalues with n = representing the largest. The example shown 
here is derived by performing PCA on 4000 simulated PSFs. As n increases, the eigenPSF tends to possess higher-frequency features. 
Interestingly, the first several eigenPSFs resemble the basis functions obtained by the shapelet formalism remarkably. However, asymmetry 
is also evident, which reveals the characteristic pattern of the data. 
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Fig. 10. — Simulated LSST image. We display only a ~ 100" X 100" cutout to show the detail. The left panel image is the result prior 
to noise addition, about 2 mag deeper than the median depth (6000 s, ~ 27.5 at 5 a) shown in the middle panel. The nominal 15 second 
exposure case is shown in the right panel. 
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Fig. 11. — Ellipticity distribution in simulated LSST images. We display "whiskers" overlaid on the LSST focal plane that we use for 
the input PSF generation. Apart from the magnitude, the PSF pattern is similar to the one in the case when only optical aberration and 
focal plane height errors are considered (Figure [81. The reduction in the magnitude demonstrates that the atmospheric turbulence mostly 
circularizes the PSF rather than introduces additional anisotropy. Because we simulate the sky near the zenith, the ellipticity due to the 
atmospheric dispersion is very small in the current case. Note that the correlation of the ellipticity change with the focal plane height 
variation is still visible. 




Fig. 12. — Residual Ellipticity distribution. Wc subtracted the model cllipticity from the observed (simulated) ellipticity. Hence, the 
length of the whiskers shows the magnitude of the residual ellipticity (5e+ , <5e x )■ 
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Fig. 14. — PSF cllipticity spatial correlation. We examine the correlation function £+(r) for PSF ellipticity. First, PSFs are paired, and 
then the tangential components with respect to the line connecting the two PSFs are evaluated. The amplitude of £_ (r) is similar to that 
of £+(r), and is omitted for clarity. Note the dramatic reduction of the spatial correlation for residual ellipticity of the stars. 
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Fig. 15. — Circularization of LSST PSF from field rotation. We illustrate a toy-model observing scenario, where a same patch of the 
sky is visited 100 times with randomized roll angles. The focal plane height fluctuation is assumed to remain constant during the entire 
mission whereas the atmospheric turbulence changes with time. Post active optics residuals are small by comparison. Wc display only the 
fraction of the focal plane for easy viewing. The observed PSF cllipticity (upper right) is much smaller than seen in individual single 15 s 
observations (upper left). The lower left panel shows the reconstructed PSF by first performing PCA on each exposure and then stacking 
all model PSFs with correct application of field rotation. The residual PSF cllipticity is displayed in the lower right panel. 
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FlG. 16. — Same as Figure [l4| but for the case that a same patch of the sky is visited multiple times with randomized roll angles. Stacking 
different PSF pat tern s reduces the correlation in intrinsic PSF (left). The residual PSF correlation function (right) remains at the similar 
level as in Figure [T4| which is not surprising because the amplitude is already at the statistical limit. 
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FlG. 17. — Null test with galaxy ellipticity correlation for a single field of view and for 20,000 sq. degrees. If anisotropic PSF effects 
are corrected, the correlation amplitude should simply reflect the statistical noise as observed. As the number of exposures increases, the 
correlation amplitude goes down because the residual PSFs become more uncorrelated (systematic noise) and also the increased depth gives 
higher number of galaxies usable for shear measurements (statistical noise). The number density of usable galaxies for weak-lensing shear 
extraction in one, 10, and 100 exposures are 10, 24, and 42 galaxies per sq. arcmin, respectively. The area of the simulated image is ~ 10 
sq. degrees. If we extrapolate the current result to the case, where all 20,000 sq. degrees of the area are taken, the expected amplitude of 
the residual galaxy shear correlation is below ~ 10~ 8 (cyan). 
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Number of goloxies 

Fig. 18. — Total (systematic + statistical) errors in galaxy shear correlation. The power-law fit gives a slope —0.97 ±0.06, consistent with 
the oc 1/n behavior of the statistical noise (solid line). A systematic noise, although not detected in the current 10 sq. degree simulation, 
will determine the minimum noise level when the nominal 20,000 sq. degree survey is completed (~ 3 billion galaxies). The dotted line 
shows the modification when we assume a 10 -8 systematic noise floor. The expected cosmological signal is greater than 10 -6 over these 
angular scales. 
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TABLE 1 

Focal plane error budget (peak-to-valley) 



Source 



Specification (Maximum Values) 



Values used for simulation 



Global focal plane adjustment 
Raft height error 1 
CCD-to-CCD height variation 
CCD tilt 

Potato chip effect 2 



10 /urn 
4.8 /jm 
4 fim 
2.5 X ltr 4 rad 
4 fim 



10 (im 
6.5 (im 
10 fim 
2.5 X 10~ 4 rad 
5 fim 



1 3 X 3 CCDs will be attached to one raft. 

2 1-2 cycles within a CCD 



